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Q^ ■ Abstract 

^^ \ We investigate the contribution that instantons make to the QCD chiral con- 

tH- \ densate for Nf = 0, 1 and 2 quark flavours. We use a simphfied model: the 

instantons are a (weighted) gas, the fermionic integrations are restricted to 
the subspace spanned by the would-be zero-modes and only the dimensional 



> 



(^ \ and chiral features of the Dirac operator are retained. Previous work along 

^^ I these lines found a power-like divergence of the Dirac spectral density at zero 






eigenvalue. By changing the various components of the model, we show that 



pi i' a power divergence appears to be a generic feature of instanton mixing for 

both Nf = and Nf ^ 0. In the latter case the divergence disappears as 
m — > in such a way that {ipip) remains finite. We find that the exponent 



X 

d I of the power decreases for larger instanton densities, becoming negligible for 

very dense instanton 'gases'. We reproduce the expected m-dependence of 
the topological susceptibility, we investigate the space-time structure of the 
eigenf unctions, as a function of the eigenvalue, and we calculate the r]' mass. 
In summary: instantons appear to provide a natural mechanism for the spon- 
taneous breaking of chiral symmetry in QCD, but they typically produce a 
power-like divergence in the Dirac spectral density. This divergence is non- 
standard but 'harmless' in the case of full QCD, and implies a pathology in 
the case of quenched QCD. 



I. INTRODUCTION 

The spontaneous breaking of chiral symmetry, and the dynamics that drives it, remains a 
major unsolved problem in QCD. However it has long been recognised |l]-0 that instantons 
provide a potential mechanism for this symmetry breaking. This can be simply motivated 
using the well-known Banks-Casher formula |Q and the Atiyah-Singer Index theorem 0. 
The Banks-Casher formula relates the value of the chiral condensate to the density of modes 
of the Dirac operator, z7(A), at zero eigenvalue: 

(^^) = 7r77(0), (1) 

while the Atiyah-Singer Index theorem tells us that the eigenvalue spectrum of the Dirac 
operator will contain 

Q[A] =n_-n+ (2) 

exact zero modes, where Q[A\ is the winding number of the gauge field configuration and 
n± are the number of zero eigenmodes with positive/negative chirality. So if we neglect 
interactions amongst (anti) instantons we get an exact zero-mode for each topological charge 
in the gauge field. This contributes a term oc 5{X) to I^(A). That is to say, the eigenvalues 
are clustered about zero. This is in contrast to a purely perturbative approach where, if 
we neglect interactions, the perturbative vacuum becomes free and the eigenvalue spectrum 
grows as A'^. That is to say, the eigenvalues are far from zero. Introducing interactions will 
shift both these distributions in ways that are not easy to predict once the interactions are 
strong. Nonetheless it is plausible that if what one wants is a non-zero density of modes 
near A = then instantons provide a much more promising starting point than the classical 
vacuum. 

An old approach to obtaining chiral symmetry breaking from instantons is as a self- 
consistent solution to semi-classical gap equations (see, for example, 0). An alternative 
approach is to perform a direct calculation of {ipip) in some semi-classical approximation 



(see, for example, [0,^0). Both approaches make drastic approximations, but the latter 
has the attraction that these can be easily made explicit. In this paper we develop upon 
the calculations which were briefly summarised in [0]. We now briefly outline our general 
approach. 

To begin, we must make more precise what it might mean for chiral symmetry breaking 
to be "due to instantons". First, it is not at all obvious that the topological charge of a 
typical gauge field must be in a form that resembles an ensemble of instantons. We will, 
however, assume that it is. That is to say, we assume that it can be decomposed into localised 
lumps of unit charge, Q = ±1, which can be characterised by the usual instanton collective 
coordinates. We shall refer to these lumps as 'instantons' even though they will normally be 
far from the minima of the action. There is some evidence from lattice simulations for this 
assumption (see, for example, P,^). A second problem is that even if the topological charge 
can be decomposed into an ensemble of 'instantons', then the correlations between these 
topological charges will be influenced by other fluctuations in the gauge fields; for example, 
those that induce confinement. So when we calculate the Dirac spectrum associated with 
the instantons (see []TU| for an example of such a calculation and |TT| for a review of related 



work) this will inevitably contain the influence of other dynamics. Thus, even if we find 
that we do have chiral symmetry breaking, it will be unclear to what extent it is due to the 
instantons per se and to what extent it is due to this other dynamics. The simplest way 
to finesse this ambiguity is to ignore all the correlations and see whether an uncorrelated 
ensemble of instantons will, by itself, drive chiral symmetry breaking. If so, if we can show 
that the spontaneous breaking of chiral symmetry is a generic consequence of the most 
general features of an instanton ensemble, then we can plausibly claim that it is indeed 'due 
to instantons'. 

Since we are interested in the implications of instantons that are generic rather than 
due to some special details of the approximations used, it makes sense to avoid details 
altogether, as far as possible, at least as a first step. Such a minimalist approach, that 
attempts to include only the essentials of the physics, was developed in 0Jl^ and we shall 



use what is almost the same model here. For the pure gauge theory (i.e. Nf = or quenched 
QCD) we shall take our instantons to have random positions and we give them some fixed 
size p. The latter embodies the anomalous breaking of scale invariance, while the former 
encodes the long distance clustering properties of the field theory. We leave the density 
of instantons as a free parameter and we ignore the other collective coordinates. If we 
think of a single topological charge, its Dirac spectrum will contain a zero eigenvalue and 
other, non-zero, modes. The former is independent of how deformed the instanton is; the 
latter depend on all the details. So if the symmetry breaking is to be a generic feature of 
the instanton ensemble, it should depend only on the zero modes. We therefore restrict 
the fermionic path integral in the background of a particular gauge field to the subspace 
spanned by the zero-modes of the instantons in that gauge field. That is to say, we are 
investigating whether chiral symmetry breaking is driven by the would-be topological zero 
modes. While the existence of a zero mode is independent of detail, its functional form 
is not. Since we do not know what this form is (and there is no reason to think that the 
classical form is at all relevant) we will try a number of different functional forms and see 
what, if anything, is insensitive to the particular choice. In a matrix representation of the 
Dirac operator, the overall structure is determined by the fact that Q = ±1 zero modes have 
opposite chirality and that 75 anti-commutes with ip. We shall assume that this is all that 
is important and simplify the calculation accordingly. We shall use the same approach for 
full QCD (i.e. Nf 7^ 0) except that we weight the instanton ensembles with the appropriate 
power of the fermionic determinant. This determinant will also be calculated in the basis of 
would-be zero-modes as described above. 

We remark that the above may be thought of as a 'bottom-up' approach to this whole 
question. There is an alternative approach that one may think of as 'top-down'. Here one 
uses random matrix theory (for recent reviews see |]T3|) identifying the appropriate univer- 



sality class on the basis of the symmetries in the problem, and then making corresponding 
predictions for the fermionic physics. To what extent the two approaches are in agreement, 
even for quantities for which they both make predictions, is not clear to us, but this is a 



question that needs to be addressed. 

The conclusion of |]7|,p^ was not only that instantons generically produce a non-zero 



chiral condensate, but that in quenched QCD this condensate diverges as m — ;► 0. This 
divergence is a reflection of a corresponding power-like divergence, 

z7(A) = a + 6/A^ (3) 

as A — i> 0, in the spectral density. The value of the exponent was found 0Jl^ to be c/ ^ 1/2. If 
this is indeed the case, then it suggests a striking pathology in quenched QCD. An essential 
first step in assessing this result is to determine how robust it is against changes in the 
assumptions employed. In particular, the zero-modes were modeled, in 0,0, by simple 
hard-sphere wave-functions (and the wavefunction overlaps were generally calculated using 
only the space-time points at which the topological charges were located). In the present 
paper we shall compare calculations with several different functional forms (and we shall use 
the full overlap functions). We shall also perform a detailed study of finite- volume effects 
and of how things depend on the instanton density or 'packing fraction' /. The latter has a 
natural definition for 'hard sphere' wavef unctions: 

where N is the mean number of topological objects in the gas, Vj is the mean volume of such 
an object and V is the volume of space-time. (Recall that the volume of a 4-sphere of radius 
p is Vi = tt'^p^/2.) What is the corresponding definition of / for other wavefunctions is not 
immediately obvious and is one of the questions we shall address below. Our calculations 
will lead us to a somewhat more nuanced conclusion than that drawn in |^ . In particular we 
find that while a divergence is indeed a characteristic product of the mixing of the would- 
be zero-modes, the value of the exponent depends quite strongly on the instanton density, 
becoming negligible for large packing fractions. 

In full QCD there should be no such pathology, and so it is of interest to ask what 
happens to this small- A divergence as we make Nf ^ 0. We find that the divergence is 



still there but that its coefficient — i> as the fermion mass — > in (apparently) such a way 
that the chiral condensate has the conventional smooth chiral limit. Thus the divergence is 
'harmless' although the Dirac spectrum is unconventional and the usual Banks-Casher Q] 
relation, given in equation |1|, need no longer be valid. 

In the next section we provide a brief discussion of instantons and chiral symmetry break- 
ing which will also serve to introduce our notation. We then introduce our minimalist model 
for an ensemble of instantons, their associated zero modes and our matrix representation 
for the Dirac operator. We point to some calculational difficulties and how we deal with 
these by the choice of an appropriate space-time manifold. We then turn to the results of 
our calculations for quenched QCD. Our main result is that the Dirac spectrum diverges 
at small eigenvalues. How this result depends on details, both unphysical - such as the 
functions chosen for the zero-modes and the space-time manifold in which they exist - and 
physical - such as the instanton density and the physical volume - is investigated in detail. 
We also investigate whether the eigenfunction of the Dirac operator changes character as the 
corresponding eigenvalue decreases. We find that indeed it does: as A — * the eigenfunc- 
tion becomes much more non-local. Having acquired some confidence that the qualitative 
features of our calculations are independent of the specific model assumptions, we then 
proceed to the (computationally) more difficult case of full QCD. Here we ask whether in- 
stantons spontaneously break chiral symmetry, what happens to the A — > divergence in the 
Dirac spectral density, and whether our model can correctly reproduce the ?Ti-dependence of 
the topological susceptibilty that is predicited by the anomalous Ward identities. We also 
calculate the rj mass and conffim that it is not a Goldstone boson. 

Some of our early results were summarised in [jl4[. A calculation of the Dirac spectrum 
using the cooled instanton ensembles obtained in the SU(3) lattice gauge theory calculation 
of p] is to be found in [|T0|. Although the results differ, in the fine detail, from what one 
obtains with a randomly placed ensemble of instantons, the main features of spontaneous 
symmetry breaking and a divergent quenched chiral condensate appear to be present. 



II. CHIRAL SYMMETRY BREAKING AND TOPOLOGY 

The usual order parameter for chiral symnietry (breaking) is the chiral condensate in the 
hmit of infinite volume and zero mass: {ilj'ip) = limm^o limy^oo {4'i')mV ^^ '^^^ express 
{ip'4')m V ^^ terms of the normalised Dirac spectral density, z7(A, m), as 

Here we have chosen to separate from V the (5-function contribution of the exact zero modes, 
whose mean number is Nz = \Q\ oc y/V. In the thermodynamic limit this zero- mode 
contribution will clearly disappear. Thus one might expect that at finite- V^ the eigenvalue 
spectrum restricted to Q = field configurations will have smaller finite- \^ corrections than 
if we allow Q to vary. (As long as the total number of topological charges is A^ ^ 1, so 
that we would have had (Q^) ^ 1 if we had allowed Q to vary, this should not represent an 
unphysical bias on the calculation.) In this expectation (to be verified later on) most of the 
calculations in this paper will be for fields with Q = 0. 

Note that if we work in a finite volume, then we would normally choose it to be at least 
large enough to accommodate the lightest particle in the theory, the Goldstone pion, i.e. 
Vi ~ 1/m^. If we insert mlf^ = m{il)il)) and crudely estimate N z — ViQ'^) = VxtV^ 
where Xt is the topological susceptibility, then we find that the zero-mode contribution in 
equation ^ is 



Ny y? — — 

_A ^ ^(^^) ^ 4(^^) (6) 

where we have used the standard values, Xt ~ (ISOMeV")^ and /,r ~ 93MeV. This demon- 
strates that the finite-volume correction to the quenched QCD chiral condensate from the 
exact zero-modes can be large even if one uses space-time volumes that might appear to be 
adequate. 

In full QCD the above calculation will differ in that the topological susceptibility will vary 
with the quark mass: Xt — ml^ipip) / ni if we have spontaneous chiral symmetry breaking. 
Inserting this into equation we find that the zero-mode contribution 



- - oc m2 (7) 



r-^ 



mV mV (m\/)5 A^/ 
vanishes when we take the chiral hmit while maintaining V large enough to accommodate 
the Goldstone pion. Note that the case of a single flavour, Nj = 1, is special. Now there 
will be no Goldstone pion because the Ua^^) symmetry is anomalous: the would-be pion is 
a massive t]'. So there is a finite mass-gap in the theory as m ^ 0. In such a situation one 
would normally expect that a fixed finite volume, with V^ > l/m,y', would be large enough. 
However we note that this would lead to a oc l/i/m divergence in the chiral condensate. 

One can ask what happens to the zero-mode contribution if one chooses to take the chiral 
limit in a fixed volume. In that case we expect symmetry restoration for sufficiently small 
m, and for smaller masses we may expect Xt oc m^f . Thus for Nf > 2 there will be no real 
m ^ divergence in {4''ip)^y at fixed V. 

The above finite volume corrections, due to the |Q| exact zero modes, induce a trivial 
divergence in the Nf = chiral condensate if we take the m ^ chiral limit at fixed 
volume. As remarked earlier, we shall evade this divergence by usually working with Q = 
field configurations. However, as described in the introduction, there is evidence 0,|T2|JT^ 
for a much less trivial divergence; one which survives the thermodynamic limit and which 
appears as a divergence in the spectral density of the form in equation |^. One of the main 
aims of this paper is to determine whether this divergence is really there, in both quenched 
and full QCD, or whether it is an artefact of the approximations used in the calculations of 
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There is a simple argument 0,0, briefly described in |^, that provides an apparently 



natural explanation for such a divergence. It goes as follows. 

In a finite volume V one expects to lose spontaneous symmetry breaking once the explicit 
symmetry breaking term is small enough i.e. once mV becomes small. That is to say, we 
expect a suppression in the finite volume spectral density for |A| < 0{1/V). Let us assume, 
for simplicity, that the only finite-V correction affecting the spectral density of the non- 
zero modes consists of a cut-off for |A| < c/V. Then the only other finite- V^ correction 

8 



concerns the exact zero-modes. If we increase the volume from V to V then their number 



increases from n^ ~ \/ {Q"^) ~ Vxt^ to ~ \/xt^- Therefore, as a fraction of the total 
number of modes, their number will decrease from ~ yXtf^ to ~ yXtP^' ■ The integrated 
spectral density (including the zero modes) does not change with volume so these two finite- 
ly corrections must match: 

/%(AMA=y|-^. (8) 

■yT 

where '^'(A) is the spectral density for V = oo. It is easy to see that this relation implies 
that 

z7(A) V A~^ (9) 

and hence that (ipip) oc m~2 in the thermodynamic, chiral limit. 

The important assumption in this argument is that as we increase the volume, the frac- 
tional decrease in the exact zero-modes is matched by an increase amongst the very smallest 
of the non-zero modes. We shall refer to this as the "eigenvalue replacement hypothesis" and 
shall investigate it in more detail later on in this paper. For now we simply point out that 
intuitively it is not unreasonable. In a volume V we expect on the average CqW exact zero 
modes. Now consider this volume split into two equal halves. Considered in isolation each 



half- volume would possess ~ coy/V/2 exact zero- modes. So in joining these two half- volumes 



together we have decreased the total number of zero modes from ~ 2 x coyVy^ = CqV^^V 
to CqVk . That means some of the original zero-modes have mixed and split from zero. But 
if V is sufficiently large, one would expect the zero modes in the two halves to be mostly 
far from each other so that the overlap will be very small, and the eigenvalues will be split 
very little from zero. This very qualitative argument should really be formulated using a 
number of subvolumes for which one can use statistical arguments. But it suffices to give 
some motivation for what we have called the eigenvalue replacement hypothesis. 

The above argument would seem to be equally applicable to full QCD, with one modi- 
fication: now Xt oc m (as long as chiral symmetry is spontaneously broken). Inserting this 
into equation | we find 
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F(A)V(!^)'. (10) 

While this is an unconventional form for the spectral density, it does not create any diver- 
gence in the chiral condensate: 

(V^V) = lim ,lim (V'^)^y 

hm / ^-^B^dX 



2{a+^} 

n — - — dy 

y2 + l 





oo 



= 7ra + 2b I -^-^ . (11) 

Jo 2/2(2/2 + 1) 

Here we have used a spectrum V{X) = a + h{m/X) 2 and have performed a change of variable 
to y = X/m. We see that the result for {ip'ip) is finite. However the naive application of the 
Banks-Casher formula, which would have produced the result an, breaks down. So, although 
a spectral density of this kind does not immediately lead to any pathological physics, it is 
clear that it might undermine any argument which relies on the A ^ and m ^ limits 
commuting with each other. 

The form we conjecture in equation |10| can also be made plausible as follows. Since a 
mode oiip will not be (individually) suppressed by the (det ip)'^f = Yli'^l'^^)^^ weighting 
if A ^ m, we might expect that the spectral density for full QCD will have the same form 
as in quenched QCD for A ^ m. If the quenched spectrum behaves as I/A2 as A — ^ then 
the spectrum of full QCD should have the form V{X,m) = {fi/X)^ for A ^ m where fi will 
have dimensions of mass, just on dimensional grounds. The theory has two dimensionful 
parameters: the quark mass m and the pure gauge scale Aqcd- Since the term (/i/A)2 
should only become important (compared to the constant term say) for A <^ m, it must be 
that fi oc m. That is to say, the form in equation |T0|. 

If we use the argument in the previous paragraph, then we can reverse our earlier argu- 
ment in order to show that in QCD we should have (Q^) oc mV. This is of interest because it 
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is actually hard to find a simple explicit argument that gives this result. (The usual deriva- 
tion using anomalous Ward identities may be simple but is hardly explicit.) For example, if 
we take the dilute gas limit then each topological charge has its separate fermionic weight- 
ing oc {detip[Ai])'^f oc m'^f that accompanies the factor of V (the charge can be located 
anywhere in the space-time volume). One immediately concludes that {Q"^) oc m^fV. Sup- 
pose instead we try the opposite extreme, where the only important effect of the fermionic 
weighting is that the gauge field acquires a factor oc m^^l'^L If we had just a pure gauge 
weighting, then we would have (Q^) oc V, and for Q^ <^V the weighting will be almost inde- 
pendent of Q. For such Q the only non-trivial weighting is our assumed fermionic weighting 
oc rn^f^'^^ = exp{— A^jIQI \n{l/m)}. It is easy to see that this implies {Q'^) oc m^f . (Now we 
do not even have the correct factor oi V.) It is clear that the real situation is much more 
subtle. When we have chiral symmetry breaking the spectral density remains finite all the 
way down to A = 0. So the \Q\ exact zero-modes may have been traded for modes that were 
arbitrarily close to zero. It is this effect that presumably reduces the naive {Q"^) oc m^fV 
dependence (which is correct when we have chiral symmetry restoration and no non-trivial 
density of modes arbitrarily close to zero) to the weaker (Q^) oc mV dependence. It will be 
interesting to see whether we can reproduce this subtle trade-off within the simple model 
that we shall be using in this paper. 

III. THE MODEL 

In this section we describe how we are going to calculate the contribution of instantons 
to the Dirac spectral density. This requires making a number of radical simplifications. 
Although we shall try to give each of these simplifications some motivation, it is important to 
perform some check that the final conclusions are robust against changes in the assumptions 
being made. 

We shall first describe the general framework of our calculation. (We refer the reader 



to []TU| for a detailed discussion of some aspects of the model.) We shall then describe the 
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various zero-mode wave functions that we use in our calculations. Comparing the results of 
these calculations will provide a check of the robustness of our scheme. 

A. The general strategy 

We are interested in determining what are the generic physical effects of instantons. 
Suppose we have a gauge field which contains rii instantons and rij anti-instantons. Suppose 
that the fields of the instantons and anti-instantons, if in isolation, would be Al and A~ 
respectively. Using a linear addition approximation one can write a field 

"7 nj 

^ = EA^ + E^7- (12) 

i J 

that represents the instanton content of the original gauge field. (We suppose that each 
individual instanton field is in a suitable singular gauge so as to be as local as possible.) 
One could do this for a realistic ensemble of (lattice) gauge fields, as was done in [10|. 



However in that case there will be correlations amongst the instantons that are partly 
determined by other, non-topological fluctuations and if we are interested in the generic 
properties of intantons per se then we do not want to include these. So the fields we 
shall consider will be those due to a randomly positioned ensemble of instantons. Such an 
ensemble clearly possesses the long-distance clustering properties of the field theory. We 
shall give the instantons a specific size; this embodies the anomalous scale-breaking in the 
theory. The instanton size sets the length scale of the theory and the volume of space-time 
can be expressed in those units. Varying the average distance between instantons amounts 
to varying the density, or packing fraction, as defined earlier. A slightly more elaborate 
approach would be to give the instantons a distribution of sizes; in particular for small sizes 
the distribution is calculable in perturbation theory. For the results of such an approach we 
refer the reader to |l^ . 



The above is plausible for quenched QCD, but breaks down for full QCD where the 
correlations between topological charges will become of increasingly long range as m — >^ 0. 
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This can be incorporated by weighting the instanton ensemble by a factor oc (deti^)^-''. 
This will be described in detail later on in the paper. 

There are a total of N[A] = rij + rij objects in our gauge field A (in the volume V), 
and the winding number of A is Q[A] = rij — rij. We know via the Atiyah-Singer Index 
Theorem [^ that the Dirac operator in this background ip[A] will have at least Q[A] exact 
zero eigenvalues. If the objects in the superposition were non-interacting then each would 
have a corresponding zero eigenvalue, 

tm^n^to = 0, (13) 



again due to the Atiyah-Singer Index Theorem. We see from equation |T^ that these zero 
modes are modes of definite chirality, in fact 7^|V'io) = ±\il>^). Interactions between the 
objects however, lead to a reduction of exact zero modes from N[A] to only Q[A]. The 
remaining N[A] — Q[A] modes have split symmetrically about zero 0J^,0. We shall fre- 
quently refer to these modes as the "would-be zero modes". The model, which is by no 
means new, is based around the idea of constructing an explicit matrix representation of the 
Dirac operator, using only the zero mode wavefunctions {|'?Aj{)), li'jo)} as a basis. Of course, 
the claim is not that these zero mode wavefunctions in isolation form a basis (they certainly 
do not span), the claim is that by constructing the Dirac operator on this subspace we can 
understand that part of the Dirac spectrum which arises from instanton interactions. Fur- 
thermore as these are the eigenvalues which generate a spectral density around zero, we may 
be able to comment upon the idea of instantons as the source of chiral symmetry breaking 
in QCD. We drop the subscript zeroes on the would be zero mode wavefunctions from now 
on, as these are the only wavefunctions which we shall consider. It is easy to see that: 

{^t\mA]\^7) = V^r (14) 

The first two identies follow from {i^,7^} = and the definite chirality of the zero mode 
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wavefunctions. The function Vij can in principle depend upon the size, position and relative 
colour orientation of all the objects in the field. We shall ignore the colour orientation of the 
objects for simplicity. And we shall assume that Vij depends only on the objects i and j. 
This is an approximation that becomes increasingly hard to justify as the instantons become 
denser. See (ITUI for a more detailed discussion. We also note that: 






(15) 



simply from the chiral properties of the zero mode wavefunctions. The function U will 
depend upon the size, position and relative colour orientation (which we ignore) of the 
objects. We now use equation ^ to construct an orthonormal set of vectors {\4'^), li^j)} 
using Gram-Schmidt orthonormalization: 



\^i) = sS-) l<t<J<nj. 
A representation of the Dirac operator (in this restricted subspace) is given by: 



nj 



ip = D 






V, 



ij 






nj 
rii 



where 



(16) 



(17) 



V = {R-yvs-^ 



(18) 



We can therefore construct this matrix representation for any given configuration of in- 
stantons providing we have been given the functions U{xf, pf, xf, pf) and V{xf, pf, xj, pj). 
Once we have the matrix it is a simple matter to calculate its eigenvalues and build up a 
spectral density using an ensemble of configurations. 
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The function U will depend on the functional form of the would-be zero modes. What 
this form is, in a realistic gauge field, is not known, apart from the fact that one expects it 
to be localised around the centre of the instanton and to have an extent that is related to 
the size of the instanton. (We remark that there is no particular reason to favour the known 
semi-classical zero- mode solution.) We assume that apart from the chirality, the spinorial 
degrees of freedom are inessential. We will therefore use wave-functions that are simple real 
valued functions of position, with an index, -|- or — , to label the chirality. We shall use 
several simple trial functional forms, as described in detail below. Only those features of 
the resulting Dirac spectrum that are independent of the wave-function form used, can be 
considered to be generic features of instanton physics. 

The function V depends not only on the form of the would-be zero modes but also on 
what we choose for ip. If all the other fields have a negligible presence in the region where 
ipf{x) and ip'i^ix) are non-zero, then we can write 



1—-^ 






iip[At + AiW^i) 
iI/)[Al]+i]/)[Ai]-i^\i^-) 
-im^) (19) 



Here the first approximation is a consequence of the linear addition ansatz as stated previ- 
ously. The second approximation holds if the objects are well separated from one another, 
however it does not hold in general (especially for the high density gases to which we some- 
times apply our model). It is however necessary if we are to progress beyond the "dilute" 
gas approximation using a model such as ours. The last two equalities hold due to the fact 
that the states are zero-modes: "liPXA^ llV'z ) = 0- Of course in any realistic situation this 
matrix element will depend in a complicated way on all the other fields. We can only hope 
to elucidate those properties that are independent of this unknown dependence. In that case 
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it makes sense to begin by simplifying this operator as radically as possible; that is to say, 
we retain only its dimensionality and its chiral properties. So, apart from anticommuting 
with 75, we assume that the Dirac operator acts like a unit operator multiplied by a scale 
factor 1/p where p is a characteristic length scale of the matrix element under consideration. 
We choose to use p = ^/pipj for the matrix element Vij (so ignoring the scale provided by 
the distance between the instanton and anti-instanton): 

VijC^^^ I d^x^~{x)^t{^) (20) 

^/PiPj J 

While extreme, we note that this choice has the correct qualitative behaviour in the dilute 

gas limit: all the eigenvalues go to zero. Of course in this limit the splitting from zero will 

depend on the detailed form of the tails of the zero- modes. In the opposite limit, where the 

density becomes large, our basis of zero-modes eventually becomes complete and it is easy 

to see that our matrix V will tend to the unit matrix divided by the size p (assuming the size 

is fixed). Note that if we had tried to do better than equation ^ and had replaced iJp\A\ 

by i ^, then the eigenvalue spectrum of V would tend to that of the free fermion spectrum 

as the basis becomes complete. This is of course no reason to panic; this limit is clearly 

unphysical and in any case all our simplifying assumptions will have certainly broken down 

by then. However it explains why in our calculations, when we vary the density, we find 

that at high densities our eigenvalue spectra appear to tend towards oc 5{\ — 1/p)- 

To obtain some idea of the errors induced by this simplification, we shall in at least one 

case below, use a more realistic version of the Dirac operator. 

B. Zero-mode wave functions 

Since our calculations are in part numerical, we can only deal with a finite number of 
instantons. That is to say, we will work in a finite space-time volume. Of course, we are 
primarily interested in the ^ ^ oo limit and so we want the finite volume corrections to be 
as small as possible. This can usually be achieved by working in a periodic box, T^*. The 
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wave-functions also then need to be periodic. The calculation of the overlaps will usually be 
too complicated to be performed analytically and will then need to be performed numerically. 
This calculation has to be performed to a very high precision, otherwise our orthogonalisation 
procedure will occasionally break down. With a large number of topological charges, accurate 
numerical calculations of overlaps can become prohibitively expensive. This means that in 
practice the wavefunctions that we can use on a (hyper)torus are quite limited. 

An alternative to working on T"^ is to calculate overlaps within W^ while locating the 
centres of the instantons within a finite space-time volume V. As V ^ oo one will obtain 
the same spectrum as one obtains by taking the volume of T"^ to infinity. The drawback is 
that finite volume corrections are likely to be much larger than for a torus. The obvious 
advantage is that it is possible to do efficient calculations with a much larger range of trial 
wavefunctions. For the sake of brevity we shall refer to such calculations as being performed 



'in M^" . 



We now list the various wave functions we shall be using in this paper. 

1. Hard sphere 
The hard sphere wavefunction is defined by: 

= otherwise, (21) 

where p- is the measure of its size and x • is its centre. This is the simplest wavefunction 
we use. In M^ the overlap is simply given by the volume of intersection between two solid 
balls in four dimensional space and this is easy to calculate analytically. 

If we are in T^ and if pj < 1/2, where / is the length of the torus, then the periodic copies 
of the hard sphere do not overlap with each other. This greatly simplifies the calculation 
of the overlap matrix element f/y-; we just calculate the overlap between the j'th object and 
the i'th object as well as with all the nearest periodic copies of the i'th object, and each of 
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these calculations can be performed analytically as though in M.'^. If pj < 1/4 things become 
even simpler since the periodic copies then make no contribution. 

Thus the calculations with ipnix) can be performed efficiently in both T^ and M^. Hence 



it is the wavefunction that was used in the calculations of ||^ and in |T^. It has some 
obvious drawbacks. In even a moderately dilute gas there will be some instantons which 
have no overlap with any anti-instantons and these will produce 'accidental' exact zero- 
modes. Moreover, one can easily show that in a dilute gas there is a trivial 1/A°'^ divergence 
in the eigenvalue spectrum as A ^ 0, that arises from the two-body eigenvalue splitting that 
results from isolated (anti)instanton pairs whose overlap tends to 0. Since a divergence as 
A ^ is one of the main issues in this paper, we will need to perform careful checks to show 
that it does not have this trivial origin. 

2. Gaussian 

Gaussian wavefunctions in M^ have tails that stretch to infinity and their overlaps can 
be calculated analytically. The Gaussian wavefunction is given by: 



1 {x-xff 



(^iv-a) = ^^^^^^p I ^j2 I ■ (22) 

where x ■ is the centre and o" • is a measure of its size. 

If one wishes to work in T^ then this wave-function must be made periodic in some way. 
We shall choose to do so as follows. First consider the following Gaussian function in M^: 



1 f {x — xf — n.iy 

^,,-,,'V--7|^exp — ^ 



G{x- xf, af, n) = ,_ ^^ exp I ^^ ) , (23) 

where xf lies in the volume V = l^ and n G Z^. Then the following provides a periodic 
wave-function: 

(x|^|.) = iV ^ G(x; xf, af, n), (24) 

where A^ is a suitable normalization constant. This is clearly a smooth function on T^, since 
we have simply added all its periodic copies to our original Gaussian. The overlap of two such 



wave-functions over T^ is simply the overlap over R^ divided by the number of T^ subvolumes 
contained in M"^. Each such overlap is a sum of terms each of which consists of the overlap 
between two simple Gaussians defined on M^. This can be calculated analytically. Hence 
our reason for choosing this particular periodic extension of the Gaussian wavefunction. 

Of course the above procedure involves infinite sums. In practice, however, a Gaussian 
falls off sufficiently fast with distance that we can often ignore all but the nearest periodic 
copies, while maintaining a very good precision. When this is so the Gaussian calculation 
on T"^ is only about 3^ times slower than the hard sphere calculation. Thus while our 
largest scale calculations on the torus will be with ipn-, we shall be able to perform realistic 
comparisons with ip^ when required. 

3. Classical 
The classical zero-mode wave-function in M^ is defined by: 

where oj, is a size parameter and x • denotes the centre of the object. In the semiclassical 
limit it would be the correct wavefunction to use, but there is no reason to think that it 
plays a special role in the full quantum vacuum. It has a slow fall-off with distance and 
for this reason we do not use it on T'^. (The natural way to proceed on T^ would be to 
determine the approximate functional form of the zero-mode numerically and to work with 
that.) 

The overlap calculation is less trivial in this case than for the hard sphere and Gaussian 
wavefunctions. Using the standard method of Feynman parameters yields the following 
integral solution: 

{i-cM>c,) ='-^[' i, f^ , . (26) 



J J J 



where s^ = {xi—Xjy. The lack of a closed form expression for this overlaps makes the calcula- 
tion far slower than when we use the hard sphere (or even the Gaussian) wavefunction. One 

19 



can indeed evaluate this integral in terms of the hypergeometric function of two variables, 
but the quickest way of evaluating this function is to calculate the above one-dimensional 
integral numerically. 

In the case of this wave-function we shall also perform some calculations that approximate 
the covariant derivative by the partial derivative, i.e. Vm — {'ip^\ ~ i ^iV^r)' I'^ther than by 
the identity, as in equation ^ There are several parameterizations of this matrix element 



for classical zero mode wavefunctions [y,|l^], and we shall use the parameterization: 



1 fi f? 

PkPi {^ + ^ iPkPi ) 
where R = |x^ ~ ^71 @- We shall label these calculations as "type II classical". This will 
provide a non-trivial test of the robustness of our results against the extreme approximation 
we have made for the covariant derivative. 

IV. QUENCHED QCD (Nf=0) 
A. Robustness 

The most striking result of the early work [0,|T^ with (essentially) the above model 
was that the mixing between the would-be instanton zero-modes produced a power-like 
divergence in the Dirac spectral density as A ^ 0. The degree of divergence that was found, 
d ~ 0.5 in equation ^, fitted in well with the simple theoretical argument given in Section 
2. The calculations were, however, only performed with the hard-sphere wave-function, i/jh, 
and for a limited range of volumes and packing fractions. Clearly one needs to test the 
robustness of this divergence against varying the details of the calculation and this is what 
we aim to do in this section. 

The first question we address is: how large are finite volume corrections when we use 
the hard-sphere wavefunction on T^ ? This is the wavefunction, see equation pl|, that was 



used in [|^. Before presenting our results we need to choose some units in which to express 
the volume. We choose our units so that the instanton size, p, is 0.2 in those units. That is 
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to say, our unit of length is 5 x p. These are the units we shall use throughout the present 
paper. We shall show the results of the finite volume study for a packing fraction of unity, 
/ = 1. This is a value intermediate between the extremes of very high and very low density 
where our approach is certain to break down (as discussed earlier). It is a natural value to 
consider in the pure gauge theory, which has only one scale. Note that for / = 1 a volume 
V = 1 will contain rij + rij = 126 topological objects. In our calculations we shall, unless 
explicitly stated otherwise, use rij = rij so that Q = 0, in the expectation (as discussed 
earlier) that this will minimise finite-V^ corrections. 

So, we have calculated the Dirac spectra for various values of V at / = 1 and for (5 = 0. 
We have fitted the small A part of the spectrum to the form in equation |^. In figure |l] we 
show how the degree of divergence, d, varies with the volume V . In figure ^ we show the 
corresponding variation of the coefficient of the divergent term as a function of the volume 
(i.e. h in equation |^). We see that both quantities appear to have finite non-zero limits 
as \^ ^ oo, and that at l^ = 1 the corrections are modest; so such a volume can be used 
to extract qualitative features of the spectrum. Although we do not show the full spectra 
we remark that finite- V corrections are small across the whole range of A and not just for 
A ~ 0. We have also performed some calculations on different space-time volumes at other 
values of / and the conclusions, of these less systematic studies are in agreement with what 
we have found at / = 1. That is to say, finite volume corrections are small and, for our 
qualitative purposes, one can safely work on a space-time volume y = 1 in our units. 

The second question we ask is: how does the degree of divergence d depend on the 
packing fraction, /, as defined in equation p* We perform our calculations on T^ using hard 
sphere wavefunctions with Q = Q and on a volume V = 1. The results are shown in figure 
^. We observe that there is a rapid decrease of d with increasing instanton density precisely 
in the range of densities that is of physical interest. Thus, while we can claim that the 
existence of a divergence is a characteristic feature of the instanton induced Dirac spectrum, 
the precise degree of divergence is not something that we can claim to predict with our 
model. We note that at very high densities the divergence effectively disappears. However, 
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since our model becomes unreliable in that limit we should be cautious about inferring that 
this is the case for the true spectrum. We also note that at very low densities the value of 
d will be determined by the detailed nature of the tail of the wave-function being used. 

The calculations described above are specified in more detail in table |. We show there 
the x^/^DF (per degree of freedom) of the best fit. We have typically fitted the eigenvalue 
range A G [0, 0.3] where the best fits have x^/^df between 1 and 2. (Note that with our 
choice of units, the maximum value of A is about 5.0). In practice one can easily reduce the 
value of x^/^DF to less than unity by reducing the fitted range to about A G [0,0.1]. The 
value of d is almost completely insensitive to which of these two ranges is chosen. In either 
case, the range of A is substantial enough that the spectrum therein is determined with a 
high accuracy so that the fact that equation |] fits well provides convincing evidence that we 
do indeed have a power-like divergence as A — i> 0. 

To emphasise the last point further we have also fitted our A — > spectra with a loga- 
rthmic divergence: 

17(A) = a + 61n(A). (28) 

Typical best x^ I^df values are shown in table |. We see that such fits are quite unacceptable 
except at the higher densities where the divergence almost disappears, (i ^ 0. In that case 
the trivial expansion A^'^ ^ 1 — (iln(A) explains the inevitable compatibility of both kinds 
of fits. We remark that a logarithmic divergence is of particular interest because of its 
appearance in quenched chiral perturbation theory [|l^ and it has been recently discussed in 
the context of Random Matrix Theory [|1^ . For a brief discussion of the possible connection 



we 



between our power-like divergence and that in quenched chiral perturbation theory [[18 
refer the reader to the concluding section of [^ and of this paper. 

We now turn to two more complex questions: first, how robust are the qualitative features 
we have described above against changes in the functional form of the would-be zero-mode 
wavefunction, and second, are these features in some way special to our particular simplifi- 
cation of the Dirac operator? 
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First, the Dirac operator. One thing we can do (and the only thing we will do) is to 
see what happens if one goes one step 'better' than we have done so far and replaces the 



covariant derivative by the partial derivative as in equation |T9|. We use the 'classical' form 
for the zero-mode, as in equation ^. We then generate two ensembles (in M^) with the same 
values of / and V. (They happen to be / ~ 1 and V c^ 1 but we have not yet discussed 
how we determine / and V for anything other than the hard-sphere wavefunction.) The 



first ensemble uses our usual simplification for the covariant derivative, as in equation |2Q . 
The second ensemble calculates the matrix element of the derivative using the approximate 
form given in equation ^. We compare the resulting Dirac spectra in figure § and table | 



(with the latter ensemble under the heading "Class. 11"). We see that while there are 
some quantitative differences, as there must be, the qualitative features - chiral symmetry 
breaking and a divergence as A ^ - are certainly both there. Indeed the values of the 
degree of divergence as listed in table | are almost identical. This provides some non-trivial 
evidence that the qualitative features we are emphasising, do not in fact arise from our 
particular simplification of the covariant derivative. 

How robust are our results against changing the wavefunction? A preliminary answer 
to this question is provided in figure ^ There we plot the spectral densities obtained with 
hard-sphere, Gaussian and classical zero-modes respectively, in a volume V = 1 and with 
a packing fraction, of / = 1 for the hard-sphere. We observe that they are indeed nearly 
identical. 

Of course, this comparison skates over some points that need elaboration. Firstly, while 
the hard-sphere and Gaussian wavefunction calculations are in T^, the calculation with the 
classical wavefunction is in M'^. Are the spectra obtained in T'^ and in R^, for a given 
wavefunction choice, similar enough for such a comparison to be useful? Of course they 
must become identical as K ^ oo; but at l^ = 1? To answer this question we show in 
figure ^ the spectral densities that one obtains with f = 1, V = 1 hard-sphere ensembles 
in T^ and M^ respectively. We observe that they are very similar. We thus conclude, from 
this and some similar comparisons, that a volume of V^ = 1 is already large enough to allow 
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comparisons on these two, somewhat different space-time manifolds. 

A second point concerns the 'packing fraction' / of wavefunctions other than the hard 
sphere. We have seen that the hard-sphere Dirac spectrum varies rapidly with /, and so this 
is likely to be the case for the other wavefunctions. So ideally one would want to plot say 
d against / and see if the results were qualitatively the same for all wavefunctions. Unfor- 
tunately there cannot be a definiton of / that is really equivalent for all the wavefunctions. 
So the best we can do is to see if one can actually find values of the size parameters a and uj 
such that the divergences obtained with Gaussian and classical wavefunctions are the same 
as the one obtains with hard spheres, for the chosen values of / and V . If that is the case 
- and if so it would represent a non-trivial result - then it will be interesting to see if the 
various ensembles possess some other properties that indicate that they are of a comparable 
density. 

We therefore proceed as follows. We consider some hard sphere configurations, in a fixed 
volume V , fixed packing fraction / consisting of objects of a single size p. To be concrete let 
us choose V = 1 (the unit four torus), p = 0.2 and rij + ni = 63 + 63 which gives / = 1. We 
extract the spectral density from these configurations and fit it using the power law form 
given in equation ^. We then generate configurations using a Gaussian wavefunction of size 
a and with the same number of topological objects as before. The latter ensures that when 
(7 is equivalent to p = 0.2, the volume will be ^ = 1. We now extract the spectral density, 
do a power law fit for A ^ and extract a value of the divergence exponent d. We repeat 
this for various values of a. If we find a value for which d is the same as for the hard-sphere 
case, then we say that this is the equivalent value. (Note that when we vary a we not only 
vary / but we also vary V. However we have seen for the hardsphere wavefunction that the 
variation oidonV is weak in comparison to its variation on /.) We repeat the above using 
the classical zero-mode. If we follow this procedure we find that a = 0.074 and uj = 0.02 are 
equivalent to p = 0.2. These values were in fact the ones used to generate the spectra shown 
in figure |^. Some properties of the various spectra are listed in table |. There we label the 
various ensembles by the 'equivalent' / and V as determined above. 
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The fact that we tuned the size parameters a and uj so that they would give the same 
divergence as the hard sphere spectrum with p = 0.2 does not mean that the comparison 
in figure ^ is in any way trivial The fact that one can successfully perform such a tuning 
already shows that power-law divergences naturally arise with all these quite different wave- 
functions. We also see in figure ^ that when we perform this tuning, the whole spectrum, 
and not just the divergent piece, is quite robust against changes in the wavefunction. 

Using the above way of setting a common scale for the various wavef unctions, we can 
extend our earlier check of the V dependence of the divergent A — i> piece of the spectrum. 
We do this for the Gaussian wavefunction in figure ^ We see that the finite-l^ corrections 
are very much as for the hard sphere. In figure ^ we show how d varies with the packing 
fraction / for the Gaussian. We see that, just as for the hard sphere, the divergence becomes 
weaker at higher densities. However the agreement is only qualitative. This means that if 
we had performed our 'tuning' to the hard-sphere spectrum at a different value of / than 
/ = 1, then we would have obtained a somewhat different equivalent value of a. This is 
inevitable since different wavefunctions are certainly not equivalent in all respects. 

It is interesting to check, in some quite independent way, whether the above 'equivalent' 
sizes do really correspond to packing fractions that are very roughly comparable. Consider, 
for example, the total overlap of an instanton would-be zero-mode with all the anti-instanton 
would-be zero modes in the configuration. This will clearly increase linearly with the packing 
fraction (since the charges are positioned at random) and will be independent of V once V 
is large enough. So its value clearly provides a measure of the density that allows a crude 
comparison between different wavefunctions. We have calculated this integrated overlap for 
the hard-sphere, Gaussian and classical ensembles used in figure |] and find values 0.496, 
0.297 and 0.495 respectively. These are indeed similar enough to be reassuring. 
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B. Realistic instanton ensembles 

One might expect that a simple way to learn what actually occurs in quenched QCD, is 
to replace the artificial ensembles that we have been using, in most of which the instantons 
have a fixed size and number and are positioned at random, with the instanton ensembles 
extracted in actual lattice simulations |^,|^ . Unfortunately such an approach possesses am- 
biguities both in principle and in practice. An example of the former is the question of how 
to separate the effect of topological from non-topological fluctuations; after all, the detailed 
featrures of the instanton ensemble will be in part due to such non-topological fluctuations. 



A practical difficulty one encounters [0 is the fact that the density of the lattice instanton 
ensembles turns out to be very sensitive to the smoothing of the gauge fields, which is nec- 
essary to remove the high-frequency fiuctuations that would otherwise obscure the (lattice) 
topological charge density. Depending on how much smoothing one performs, one can obtain 
p!0| spectral densities with or without a small-A peak. 



That is not to say that calculations with realistic instanton ensembles don't have in- 



teresting lessons to teach us. For example one learns [|T0| that within an apparently dense 
instanton 'gas', the narrower instantons can contribute a small-A peak to the Dirac spectral 
density, just as though they were part of a much less dense gas. The reason for this phe- 
nomenon is simple [|I^. The high packing fraction of the lattice instanton 'gas' is dominated 
by the very large instantons (since V} oc p^ in equation ^ but such instantons will have a 
small part of their wavefunction in the small region of space-time occupied by the wavefunc- 
tion of a small instanton. Thus the matrix elements of Up between large and small instantons 
can be neglected to a first approximation, and the latter effectively interact as though they 
belonged to a much less dense gas of, on the average, smaller instantons. This observation is 
of particular interest because it (potentially) provides a justification for the use of relatively 
dilute gases of small instantons in various successful phenomenological studies |p|. 



In ||T0[ the would-be zero mode wavefunctions used were hard spheres. However we have 



seen in this paper that while we can set up an approximate equivalence between different 
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wavefunctions by matching spectral densities or, more directly, the overlaps themselves, the 
actual scales characterising the matched wavefunctions are quite different. Since we have not 
determined which scale corresponds to a particular instanton size as determined in a lattice 
calculation, this again leaves a great deal of ambiguity in determining the corresponding 
packing fraction, /, since doing so would require determining the matched hard-sphere 
radius to be used in equation Ij. This is an ambiguity of which the reader of flOl should be 
aware. 

Indeed when we refer in this paper to the size of an instanton, for example p = 0.2, the 
reader should be aware that what we mean is that the instanton has such a size that its 
would-be zero mode can be approximated by a hard-sphere of radius p = 0.2. (That the 
real would-be zero-mode can be approximated by a hard sphere is suggested by the fact that 
we can match the spectral density from hard spheres with that from Gaussians or classical 
zero modes.) Its actual size, in terms of the spread of the topological charge density, may 
well be somewhat different, and trying to determine that (for the fully fluctuating vacuum 
of interest) lies well beyond the scope of the present work. 

C. Dynamics 

We have provided evidence that the power-like divergence induced in the Dirac spectrum 
by instantons is not an artefact of the approximations used in our model but is a real effect. 
We can therefore begin to ask some questions about its dynamical origins. There are three 
issues we shall address. The first is whether there is any evidence for the simple picture for 
the origin of this divergence that we presented in Section |I|. The second question is whether 
the divergence has an origin in the simplest two-body interactions between instantons and 
anti-instantons. Finally we ask whether the eigenvectors of the Dirac operator have any 
striking features as A ^ 0. 
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1. Eigenvalue replacement hypothesis 

In Section ^ we described a simple scenario in which a divergent spectrum will naturally 
arise. In that argument one assumes that the only finite-^ correction in the (normalised) 
Dirac spectrum consists of a cut-off for |A| < c/V (for some constant c). If we increase the 
volume from l^ to V + 5V then the (fractionally) greater number of modes in the interval 
|A| G [c/(l^ + 5V),c/V] must match the fractional decrease in the number of exact zero- 
modes. This assumption that the fractional decrease in the exact zero modes is compensated 
for by an increase amongst the very smallest non-zero modes, we earlier referred to as the 
'eigenvalue replacement hyopothesis'. As we saw, it leads to a divergent spectrum, as A — i> 0, 
with a power exponent (i = 0.5. The fact that in practice we find that the power depends 
on the packing fraction tells us that this argument is approximate at best. It is interesting 
nonetheless to see if it does capture some of the essential features of what is happening. 

A simple first step (which is as far as we shall try to go in this paper) is to determine 
whether our eigenvalue replacement hyopothesis is qualitatively correct. Since this concerns 
the nature of the finite-V^ corrections, we shall need to generate configurations with no 
constraint on the value of Q. Our strategy will be twofold. Firstly we generate configurations 
at a fixed packing fraction but for two volumes, and we check explicitly whether the finite 
volume corrections to the spectrum are really concentrated at small values of A. Secondly 
we shall perform calculations at fixed V and compare, for example, configurations with 
(5 = and Q = 2 so as to see whether the extra 2 non-zero modes in the former are indeed 
concentrated amongst the lowest modes of the latter. 

For our first calculation, the one with two different volumes, we draw the number of 
instantons and anti-instantons in each configuration from a normal distribution: 

Ni^A ~ N{N/2, iV/4) (29) 

where iV oc V^ is the mean number of topological charges in the volume V (for instance 
we have N = 126 for V = 1, / = 1). We see therefore that the total number of objects 
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N = Nj + Na ~ N{N, N/2) and that the winding number distribution follows Q = Nj — 
Na ~ A^(0,iV/2). (The actual numerical factors arise as this is nothing other than the 
central limit theorem applied to the binomial distribution Ni^a ~ B{N, 1/2) for large N.) 
The resulting Dirac spectra for such a gas in two different volumes is shown in figure |^ and 
is listed in table |I|. The results clearly show a far greater difference at small eigenvalues 
with the curves converging for A > 0.04. (Note that the spectrum extends out to Xmax ~ 5.) 
This can be seen more clearly if we simply focus on the difference between the two spectral 
densities as in figure H. This provides qualitative support for what we have called the 
eigenvalue replacement hypothesis. 

As an aside, we note from table H that the value of d for ensembles where Q varies 
shows much stronger finite-^ corrections than for the Q = ensembles listed in table |. (As 
V — i> oo both ensembles should give the same value of d. ) This shows that our earlier strategy 
of using Q = ensembles so as to minimise finite volume corrections was appropriate. 

One way to go beyond the above qualitative analysis is as follows. Consider configurations 
with some fixed total number of charges i.e. in a fixed volume. We denote the un-normalised 
density of non-zero eigenvalues obtained from configurations with topological charge ±Q by 
p(A, \Q\). Note that we explicitly exclude the exact zero-modes from this density. Its integral 
over the interval [Ai, A2] is simply the average number of eigenvalues in that interval for such 
configurations. Let us now define a second density, Pmini^, \Q\)^ that counts only the smallest 
positive eigenvalue in each configuration. So its integral over A G [0, 00] is just 1. Suppose 
we now simultaneously consider configurations with charge \Q\ and \Q\ +2. The latter will 
have 2 extra zero-modes, i.e. the former will have one extra positive mode (and of course a 
corresponding negative mode as well). In this context the eigenvalue replacement hypothesis 
would say that the two extra zero- modes in the configurations with charge |(5| + 2 replace 
the lowest positive (and largest negative) eigenvalues in the configurations with charge \Q\, 
i.e. 

PrmniX, \Q\) = P(A, \Q\) " p(A, \Q\ + 2). (30) 
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We shall test the following integrated form of this relation. The smallest eigenvalue is 
obviously localised at small A. Define Aq by 

Ao 

/"p™„(A,|g|) = 0.98. (31) 



Thus the interval [0, Aq] is essentially where the lowest positive eigenvalue is localised. (We 
have used 0.98 rather than 1 in equation RT] because formally at least there will be a non- 



zero, even if infinitesimal, probability for the smallest eigenvalue to be large.) We can now 
calculate the quantity 

Ao 

5p{\Ql \Q\ + 2) = y"{p(A, \Q\) - p(A, \Q\ + 2)}. (32) 



If the eigenvalue replacement hypothesis held exactly then we would expect this to be 0.98 
for all \Q\. This test was performed in |l^ using hard spheres and we quote the results of that 



calculation in table |ITl| . (In [0,|l^] the overlap calculation is cruder than ours, but that should 
not make a significant difference.) We observe from the table that the eigenvalue replacement 
hyothesis does indeed provide an approximately valid description. Typically about 80% of 
the difference between the \Q\ and |<5| + 2 spectra resides in the region occupied by the very 
lowest eigenvalue in the configurations with \Q\. We also note that as \Q\ increases this 
percentage decreases. It is presumably this feature that leads to the power of the divergence 
not being precisely 0.5 and is what enables it to vary with the packing fraction. Thus, while 
the eigenvalue spectra do not satisfy the eigenvalue replacement hypothesis exactly, they do 
so approximately and it presumably provides a qualitative starting point for understanding 
the origins of the divergence. 

2. Two body interactions 

In the case of hard-sphere wavef unctions, when the gas is sufficiently dilute that the 
chance for an instanton to overlap with more than one other charge may be neglected, it is 
easy to show analytically that the eigenvalue spectrum diverges as A —> as 1/A°-^. This 
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follows from the fact that if two hard spheres overlap by a small distance S then the overlap 
volume is oc ^2. This will then be the eigenvalue shift from due to the two body mixing. 
The probability measure is oc d6, which, because A oc ^2, is just oc 1/A°'^. Thus here the 
origin of the divergence is relatively trivial. The fact that our observed divergence appears 
for all the quite different wavefunctions that we have used, already tells us that its origin 
can hardly be so trivial. Nonetheless it does raise the question whether it is primarily due to 
2-body dipole-like mixing, in which case it should be easy to understand, or to much more 
complex multi-body mixings. This is the question we shall now address. 

We address this issue in two ways. Firstly we generate a 'background spectral density' 
consisting of simple pairwise splittings of eigenvalues. Secondly, we explicitly construct a 
gas of dipoles and see if we still find a peak at small eigenvalues. 

Figure ^ shows the spectral density for the hard sphere wavefunction coming from the 
actual eigenvalues of our Dirac matrix. It also shows a 'background curve' which originates 
from pairwise interactions between the topological charges in the same ensemble calculated 
as follows. Consider a configuration with Q = 0. We go through the set of instantons and 
for each one find the largest overlap with a neighbouring ant i- instant on. We then use this 
overlap in a 2-body mixing matrix to obtain the eigenvalue shift from zero: 

Xi = ±m.ax{Vij} (33) 

j 

This then produces for us the 'background spectral density' plotted in figure ^ We remark 
that if we had Q < then we would again go through the rij instantons, associating each 
with two eigenvalues, symmetrically split about zero, as above. The rest of the eigenvalues 
are zero by the the Atiyah-Singer theorem. Analogously for Q > 0. The background curve 
thus generated only takes into account pairwise splitting and furthermore ignores totally the 
effects of other objects of the same charge in the vicinity. We see that the A — > peak does 
not appear in the background curve. This indicates that the A ^ peak that we observe in 
the Dirac spectrum is not due a simple 2-body interaction. 

There are many ways to produce such a background curve, so this example cannot be 
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considered conclusive. For example one might try to produce a background curve just as 
above except that any anti-instanton that is used once in a 2-body mixing cannot be used 
again; it is effectively removed. This produces a quite different spectrum, which has a 
plethora of extra near-zero and exactly zero eigenvalues. This is because, as we remove 
anti-instantons from the configuration, we are soon left with a very dilute gas of topological 
charges. This is clearly an artefact of the procedure, but serves to illustrate the inherent 
ambiguity in trying to generate a background curve. 

To further explore this question we explicitly construct a Q = ensemble of topological 
charges from a gas of dipoles, pairing off each object with one of the opposite chirality. 
We use hard sphere wavefunctions of size p and the distance between the opposite charges 
within each dipole is chosen at random (with a four- volume weighting) within some interval 
|x^ — x^l G [0, Tfnax] for somc Tmax- If fmax < 2p then the dipoles will always have a non-zero 
overlap and we would naively expect no eigenvalues near zero if only 2-body interactions are 
important. In figure |10| we show the spectrum for such a gas of dipoles in the case where 
fmax = 2p. We compare it to the spectrum one obtains with a Q = gas of (anti)instantons 
placed at random in the volume. We also show the background spectrum that one obtains 
by calculating the eigenvalues of the dipole gas, from the mixing within each dipole taken 
in isolation. (Note that this may differ somewhat from our previously defined background 
spectrum.) We observe that the dipole gas spectrum is far from the background curve and 
not so different to that of the randomly placed charges. This shows that the eigenfunctions 
corresponding to the small modes in the divergent peak, are formed by a mixing that is very 
far from simple two body. 

An even more striking example is provided by choosing Vmax = P- Now the minimum 
overlap is very large and we would naively expect the spectrum to have no small eigenvalues 
at all. In figure ^ we show the background curve obtained by treating each dipole in isola- 
tion, and indeed it possesses a sharp cut off as expected. Remarkably enough, however, the 
actual Dirac spectrum we calulate from this gas of heavily overlapping dipoles does possess 
a non-zero density of modes at A = and hence breaks chiral symmetry spontaneously. This 
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shows that there is a component to the mixing that is far from simple two-body, and it is 
this component that drives the spontaneous breaking of the chiral symmetry. 

3. The eigenf unctions for A ^ 

It would obviously be interesting to learn something about those eigenfunctions of the 
Dirac operator that correspond to the small eigenvalues that drive chiral symmetry breaking. 
We have seen above that their origin lies in something more complex than simple two- 
body mixing. Since chiral symmetry breaking is associated with the presence of a massless 
Goldstone boson, the pion, one might expect that the eigenfunctions would become more 
extended as A ^ 0; perhaps, even, that the extent of these eigenfunctions might diverge. 
If this were so, then this could perhaps be linked to the masslessness of the pion using the 
standard decomposition of the component quark propagators in terms of the eigenvalues 
and eigenfunctions of the Dirac operator ||19|| . 

To explore this question we need a measure of the dispersion of the eigenfunction \ex) 
corresponding to the eigenvalue A of our Dirac operator, i.e. i^\e\) = X\ex). We first define 
the function 

R{x,) = {ex\{x-x,f\ex), (34) 

that is clearly a measure of the dispersion around the point Xc- The extent space-time of the 
eigenfunction is provided by the dispersion about the 'centre-of-mass'. Since the dispersion 
is minimised when we choose Xc to be the 'centre-of-mass', we can calculate the extent of 
the eigenvector directly by 

D = mmR^^\xc). (35) 

In practice we do these calculations by decomposing the eigenvector Ica) into a linear com- 
bination of the would-be zero-modes and working out the dispersions of these using, where 
appropriate, some simple approximations. 
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We have performed a calculation with an ensemble of hard-spheres on volumes V = 
0.41, 1.0, 1.5. The results of -D(A) for these three volumes is shown in figure |1^. We see that 
the dispersion increases with decreasing eigenvalue; perhaps roughly linearly. We also see 
that there are stronger finite volume effects at smaller values of A. It would appear that the 
eigenvectors for small A are being constrained by our finite volumes to be smaller than they 
would otherwise be. Whether the dispersion would actually diverge at A = in the V — i> cxo 
limit is a question we obviously cannot answer. Certainly this shows that the modes that 
are important for the spontaneous breaking of chiral symmetry are delocalised on the scale 



of instanton nearest-neighbour interactions. The implications of figure |T2| are intriguing and 
need further investigation. 

V. FULL QCD (Nf=1,2) 

In this section we shall use the simplest hard-sphere version of our model to analyse 
QCD with 1 and 2 quark flavours. The calculations are, computationally, far more intensive 
than for the pure gauge theory, and we will not be able to perform the various checks that 
we were able to perform in that case. Thus, although we are reassured by the fact that this 
simplest hard-sphere analysis proved remarkably robust in quenched QCD, we must regard 
the calculations here as being of an exploratory nature. 

There is a surprisingly large number of non-trivial questions we can ask within the 
framework of our model. We shall focus upon: 

• Spectral density. How does the spectral density behave with dynamical quarks ? Do 
we still get a power divergence as seen previously, b{m)\~'^^"^' where b{m), d{m) are 
now dependent upon the quark mass ? What is the behaviour of the spectral density 
as a function of the number of quark flavours Nf ? 

• Chiral condensate. What is the behaviour of {tpilj){m) ? In particular, do we have 
spontaneous symmetry breaking and how does it behave as a function of the number 
of quark flavours ? 
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• Topological susceptibility. General arguments give the behaviour of this quantity as: 

(Q^) oc mV symmetry broken phase 

oc m^iV symmetric phase. 

Can our model reproduce such behaviour ? 

• Particle masses. In particular, does the mass of the r] meson go to a non-zero limit as 
the quark mass m ^ ? 

We shall begin by describing in detail how we generate the instanton configurations for 
Nf 7^ 0. We then describe how we calculate the r] and a masses in our model. From there 
we move on to our calculations, in which we address the questions listed above. 

A. Ensemble generation 

The weighting of the instanton ensembles will have a gauge part and a fermionic part. 
The gauge part will be the same as in the quenched case, except that we will allow both 
Q and A^ (the total number of topological charges) to vary. We vary Q because we want 
to calculate the m-dependence of (Q^). In addition we expect that {Q"^) will cease to be 
large once m is small, and then there is the danger that using only Q = Q will bias the 
physics. The reason for varying A^ is less compelling and, in practice, the main purpose of 
such calculations will be as a check on our more extensive fixed- A^ calculations. 

1. The gauge weighting 

As in the quenched case, the topological charges are positioned at random in space-time, 
with a fixed size (p=0.2). When we wish to allow A^ = A^^ + Nj to vary, we choose the 
number of instantons using a Poisson distribution, and the same for anti-instantons. The 
two distributions are taken to be independent: 

P{Na = s,Nj = t)= exp(-2/i)^ (37) 
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Note that fi would be the mean number of instantons (as well as of anti-instantons) if we 



used the gauge weighting alone. (This is similar to equation |2^ with /x = N/2, albeit with 
a different variance.) 

2. The fermion weighting 

In the case when A^ = A^^ + Nj is fixed, the obvious fermion weighting to use is one that 
is calculated entirely within the basis of would-be zero modes, 

mm{Ni,NA) 

det(z^[A]-zm)=ml^^-^^l JJ (A,' + m') (38) 

4 = 1 

(raised to the power Nf) where the Aj are the eigenvalues we obtain, for the given instanton 
ensemble, as described earlier in this paper. There are of course other modes but their 
number does not change, since the total number of modes N^ is fixed if we are in a fixed 
volume (with a fixed ultraviolet cut-off). We have assumed, in our calculations of the chiral 
condensate, that it is a good approximation to neglect the variation of these other modes as 



we alter the instanton gas, and equation |3§ embodies the same assumption. 

If, however, we vary A^ = N^ + Nj then we can no longer ignore these other modes: 
their number, Nt — N, will vary and so some choice has to be made for their weighting. We 
will make the simplest choice: each mode will be set to a fixed value Xnz- In that case our 
weighting becomes 

inin{Ni,NA) 



detmA] - tm) = {%z + m2)(^--^^-^^)/2ml^^-^^l J] (A^ + m") (39) 



/— 2 

i=l 

(again raised to the power Nf). It is clear that we do not actually have to choose a concrete 
value for A^^ since it is only the ratio of the determinants that enters into the probability 
of changing one configuration of instantons for another one, and in the ratio the factor 
i^NZ + "^^)^^^^ drops out. However we do need some choice for the value of Xnz- In- 
keeping with our assumption that the lowest-lying eigenvalues are generated by the mixing 
of the would-be zero modes, we shall choose Xnz to be somewhat larger than the average 
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eigenvalue in the zero-mode basis. A typical quenched QCD eigenvalue spectrum is shown in 



figure |T^, and this has an average eigenvalue of A ^ 1.1; so we shall choose \nz = 2.0 in our 
calculations. The somewhat arbitrary nature of this choice is unsatisfactory and in practice 
we shall minimise the ambiguity by only performing such calculations for m <^ Xnz- (We 
remark that in our Monte Carlo we have actually used {Xnz + ttT') rather than {Xjy^ + m^Y^'^; 
but for m <^ Xnz this should not matter.) Most of our calculations will be for fixed A^. In 
this case the values of Q are necessarily all even or all odd. This should not introduce a 
significant bias as long as (Q^) is not too small. 

3. Monte Carlo simulation 

As described above, we incorporate a gauge and fermion weighting into our model using 
the two parameters ^ and Xnz- The parameter /x is related to the packing fraction of our 
configurations but does not determine it wholly because the fermion determinant will also 
play a part in finding the equilibrium number of objects in the gas. 

The Monte Carlo simulation begins with a random configuration. We then move a single 
object to generate a new trial configuration. This process is repeated as we "sweep" through 
the gas, moving each object in turn, accepting moves according to the standard Metropolis 
algorithm: 

P(accept) = 1 det(new) > det(old) 

= {^iJt^}'^' otherwise. 
Periodically (usually every 10 attempted moves) we attempt to either increase or decrease 
the number of instantons or anti-instantons by one. We expect the system to come into 
equilibrium after some number of sweeps. As the change between successive configurations 
is small (differing only in the position of a single object or in having one extra or one 
fewer object), we clearly require long sequences of configurations, both to thermalise and 
to obtain a useful number of effectively independent configurations. Whilst we use all the 
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configurations that we generate (after equilibriation) our statistical analysis uses a binned 
jack-knife procedure that should ensure reasonably accurate error analyses. 

To ensure thermalisation we monitor the values taken hj N = Ni + Na and Q = Nj — Na 
during the sequence of instanton configurations generated by the Monte Carlo. In practice 
we find that a sequence of 630000 configurations appears to be sufficient to explore the 
available phase space; and that after 63000 configurations we appear to have thermalised 



(i.e. lost all memory of the random starting configuration). This is illustrated in figure ^ 
and figure |15| for the lightest quark mass we use, and in figure |16| and figure |l^ for an 
intermediate quark mass. As one would expect the phase space is explored more weakly for 
the lighter mass. 

B. Calculating masses 

The usual way to calculate masses in a (Euclidean) field theory is from the t-dependence 
of appropriate correlation functions. The latter can be estimated by standard Monte carlo 
techniques. We can attempt to apply the same approach to our instanton ensembles. Since 
we calculate the eigenvalues (and can easily calculate the corresponding eigenvectors) of 
ip, we can, in principle, calculate quark and hence hadron propagators. Here, however, we 
shall limit ourselves to calculations that involve only 'gluonic' operators. That is to say, 
we calculate the masses of fiavour-singlet mesons. This includes what is perhaps the most 
interesting meson in this context, the rj'. We shall also calculate the mass of the scalar a. 

The rj' has quantum numbers 0^^ which are also the quantum numbers of the topological 
charge density. Thus it can be calculated from the correlation function, {Q{0)Q(t)), assuming 
the usual decomposition 

(Q(O)Q(t)) = 5^c„exp(-M„t). (41) 

n 

The lowest mass should be the rj'. In practice what we do is to divide up our volume into 
a number of "strips" each of width 6t <ti I. We calculate the total 'charge' in a given strip 
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simply by adding up the charges of all the objects whose centres are contained within the 
strip. We call this Q{t) for the strip [t,t + 6t). That is to say, Q{t) = Nj{t) — N^it) where 
Ni/A{t) ai'e the number of instantons/anti-instantons whose centres are contained within the 
strip. This differs from the actual topological charge density within a strip, but is clearly an 
equally good operator. (The difference is due to the spread of the instanton core; but that 
is modelled by some relatively arbitrary function within our model and so there is no point 
in including that.) 

One might worry that the lightest mass contributing to the correlator in equation ^ 
might be a ^ glueball rather than the rj'. However we note that in the pure gauge theory 
the topological charges are completely uncorrelated, so the correlation length is zero, i.e. 
the 0""*" glueball mass is cx). Thus it seems safe to assume that it can be ignored in full QCD 
as well. 

The a has quantum numbers O"*"^ and so can be calculated from correlators of the glu- 
onic action density, {Sg{0)Sg{t)). Within our model the gluonic action density is sim- 
ply proportional to the density of topological charges and so we use the correlator of 
N{t) = Nj(t) + NA(t), i.e. the total number of charges whose centres are located within 
the strip labelled by t. For the same reason as in the case of the O""*", we can ignore the 
glueballs. However there is now a non-trivial vacuum state contribution, and so we have to 
use the vacuum subtracted correlator, {N{0)N(t)) — {N{0))^. The large-t exponential fall-off 
of this correlator should then give us the mass of the a meson. (Note that if we fix A^ in our 
calculation, this induces long-range correlations that are not due to particle exchange. For 
this reason we do not calculate the a mass in our fixed- A^ calculations.) 

To make sure that one has identified the large-t exponential decay of a correlation func- 
tion, it is useful to define an effective mass 

ejtK J V(c(o)c(t-i))y ^ ^ 

where we have written our correlator with a generic operator O which can represent either 
the charge or the number density. It should be apparent that if the correlation function 
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is in fact given by a simple exponential for t > to then our effective mass meff{t) will be 

independent of t, within errors, for t > to- If tmin is the smallest value of to for which this 

is the case, then meff{tmin) provides us with an estimate of the lightest mass, and its error. 

It is of course the case that our model is not a field theory and a decomposition of the 



kind given in equation O may not hold. However, since our model is intended to provide an 
approximation to the relevant infrared physics of QCD, we will take the masses we calculate 
in this way as being analogues, within our model, of the corresponding masses in QCD. 

C. Calculations 

We perform calculations for both 1 and 2 quark flavours. For Nf = 1 one expects no 
Goldstone boson since the one would-be Goldstone boson, the rj', acquires a mass through 
the anomaly. (For simplicity we shall refer to the flavour singlet pseudoscalar as the r]' even 
though that is, strictly speaking, the flavour singlet member of the nonet with Nf = 3.) 
Although the U{1) chiral symmetry is broken by the anomaly, it is still interesting to see 
whether the Dirac spectrum possesses a non-zero A ^ limit. The Nf = 2 case, on the 
other hand, contains essentially all the features of the physically relevant Nf = 2 + 1 case. 

We will carry out calculations separately for A^ flxed and for A^ being allowed to vary. 
The latter calculations will show that (A^) varies weakly with the quark mass m; in contrast 
to (Q^). This suggests that constraining A^ to be flxed is not a serious bias. Most of our 
calculations will in fact be for A^ flxed, with the variable A^ calculations serving as a check. 
Fixing A^ has the advantage that we can avoid our rather ad hoc representation of the non- 



zero modes through an average eigenvalue, Xnz- In addition we can see from flgures |T^ ^ 
that the Monte Carlo explores the fluctuations in A^ much more slowly than the fluctuations 
in Q. Thus a flxed-A^ calculation may also prove to be statistically more accurate. 

All the calculations described here will be with hard-sphere wavefunctions with radius 
p = 0.2 and in a space-time volume V = 1. Our flxed A^ calculations will contain A^ = 126 
which corresponds to a packing fraction / = 1.0. When we vary A^ we choose the gauge 
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weighting so that the average value of A^ would be 126 if there were no fermionic weighting. In 
practice the latter reduces (A^) so that the packing fraction lies in the interval / G [0.35, 0.55]. 

For the calculations at fixed A^, we choose quark masses m = 0.15, 0.3, 0.5, 1.0, 2.0 
and 3.0. For the calculations with varying A^ we only use the lowest 3 values of m. This 
is because the fermionic weighting will presumably become sensitive to the precise value of 
\nz once m becomes comparable to it. 

We expect that in our finite volume, we will lose chiral symmetry breaking once m 
becomes too small. To gain some intuition as to when that might occur, it is useful to 
express everything in physical 'MeV units. If we take the instanton radius to be about 
0.5/m P] then the length of our space-time, being 5 x p, is about 2.5/m ^ l/80MeV. So 
m = 1 corresponds to m = 80MeV, i.e. just a little lighter than the strange quark mass, 
and m = 0.15 corresponds to m ~ 12MeV, i.e. close to the physical u, d masses. Thus 
we should certainly not be surprised to see finite- 1^ effects at the lower end of our range of 
masses. 

We now turn to the detailed results of our calculations, the parameters of which are 
summarised in table |IV|. 

1. Nf = 1; variable N . 



In figure |I8| we display the spectral density of the Dirac operator for quark masses, 
m = 0.15, 0.3, 0.5. The first thing we note is that the spectral density appears to diverge 
as A — > for all three masses. It is however, equally clear that the spectral density is 
not independent of the quark mass. Indeed, if we fit the spectra to the power-law form in 
equation ^ then we find, as shown in table 0, that the coefficient of the divergence appears 
to vanish, 6 — * 0, as m — > 0. Thus the divergence in the spectral density need not lead 
to a divergent quark condensate of the kind that we encountered in quenched QCD. This 
appears to be confirmed by a direct calculation of the chiral condensate against m, as shown 
in figure O, where we see that we appear to have chiral symmetry breaking with a finite 
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condensate. 

We also calculate (Q^) as a function of m, and list the values in table 0. As we see 
from figure EDI the variation of {Q"^) is approximately linear with m in accord with what we 



expect from the anomalous Ward identities for QCD: 



V N] 



+ 0{m'). (43) 



Indeed, if we insert into this equation the values relevant to our simulations, i.e. V = 1, 



Nf = 1 and the value for {ipip) shown in figure |T9|, then we obtain predicted values for {Q ) 



that are within ~25% of the values shown in figure 20. This is quite remarkable 



We now turn to our calculation of the flavour singlet meson masses. In figure 21 we show 



an example of the correlation functions from which we extract the t] and a masses. We can 



use equation ^ to extract effective masses, as shown for the r] in figure ^. The errors are 
large and there is not really much evidence that meff{t) plateaus at large t to a finite mass, 
m^' . In this particular case we choose the effective mass at t = 0.1 as our estimate of m^' 
since, roughly speaking, effective masses at larger values of t are consistent with it within 
errors. We list in table [V] our mass estimates, obtained in this way. Despite the crudity of 
our calculation it seems clear that m > does not vanish as m ^ : the f] is not a Goldstone 
boson. In our units, where p = 0.2, we find m ' ~ 12. If we transform to 'MeV units, 
using p ~ 0.5fm, then this mass becomes mi ~ IGeV. In QCD this mass is generated 
by topological fluctuations, and through large- A^c arguments can be related to the value of 
{Q"^) in the quenched theory. With our gauge weighting we would obtain in the quenched 
theory {Q^) = N = 12Q and hence, in 'MeV units, {Q^)/V ~ {260MeV)^. This is not so 
far from the true value PJ in quenched QCD, ~ {200MeV)'^, and shows that our model is 
doing remarkably well in reproducing the link between topology and physics in full QCD. 

In table [V] we also list our estimates of the mass of the a. It is comparable in mass to 
the r] , although perhaps slightly lighter. 
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2. Nf = 1; fixed N. 

The results of our calculations at fixed A^ are summarised in tables ^ and 0. Since the 
packing fraction is larger than for the variable- iV calculations (126 topological charges per 
unit volume versus ~70) we do not expect agreement even in the common m-range. However 
the qualitative agreement is quite satisfactory. In figure ^ we plot the spectral density in 
this mass range and we see the same features as we saw in figure |TH|. The chiral condensate 
is plotted in figure ^. We see that there is a finite condensate as m —* 0, but that it is 
larger than the one shown in figure 110. However the factor by which it is larger is close to 



the trivial factor of ~ 126/70 (the ratio of the number of charges). As for {Q"^), we see an 
approximate linear decrease for m < 1.0. The actual value is larger than for the variable- A^ 
case although by less than the factor of ~ 126/70. (One would only expect precisely this 
factor in the dilute gas limit.) For larger values of m the linear increase must flatten off, as 
it does, because it must approach the quenched (m — > oo) value of {Q^) = N = 126. The 



calculated values of the r] mass, as shown in figure ^, are consistent with those obtained 
previously (and confirm that the rj is indeed not a Goldstone boson). All this confirms our 
expectation that fixing A^ does not alter the physics of the model in any significant way. 

The larger range of m gives us a clearer perspective on the oc b/X'^ divergence in the 
spectral density. We see from table |V| that the power d is independent of m except at the 
very smallest value of m. We provisionally ascribe this to a finite volume effect. (The reason 
for doing so will become clearer once we present the Nf = 2 calculations.) We note that 
the value of d is smaller than the one we obtained for the variable- A^ calculation. Since the 
packing fraction of the latter is smaller, this mirrors what we found in quenched QCD: the 
divergence weakens with increasing instanton density. Although d appears to be independent 
of m, this is clearly not the case for the coefficient b. In fact, in the intermediate mass range, 
m E [0.3, 1.0], the dependence of b is consistent, within errors, with b{m) / X'^^'^^ = bQ{m/XY° 



as suggested by the simple argument leading to eqn 10 
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3. Nf = 2; variable N. 

The spontaneous breaking of chiral symmetry that we found for Nf = 1 is academic 
since the symmetry is in any case anomalous. We now turn to the corresponding Nf = 
2 calculations where there is a non-anomalous part of the chiral symmetry that may be 
spontaneously broken. 

Before proceeding, a note of warning. It is of course trivial to vary the number of fermion 
flavours: all that is required is to take the ratio of the determinants to the power Nf in the 
Metropolis step. Not surprisingly, however, this makes the exploration of the available 



phase space much slower. This is illustrated in figure ^ where we show the total number 
of topological objects in a Monte Carlo generated sequence of configurations at our lowest 
mass value, m = 0.15. The corresponding plot for Nf = 1 was shown in figure |l^. The 



variation of A^ is very much slower in the case of 2 flavours. It is clear that if we want to go to 
smaller masses or larger values of Nf we shall need either much longer Monte Carlo sequences 
or improved algorithms. This point is further emphasised by flgure ^ which displays the 
increasing difficulty in moving between sectors of different net topological charge. 

The Nf = 2 calculations with variable A^ have been performed for masses m = 
0.15, 0.3, 0.5 (chosen, as before, so that m <^ Xnz)- The spectral densities for these 
three masses are depicted in flgure 123. A quick comparison with flgure ITS indicates that 



something very different is occuring for two flavours of fermions. The spectral density is 
much smaller in magnitude, and decreasing rapidly as m ^ 0: all the signs of chiral symme- 
try restoration. This is highlighted by a direct comparison of the Nf = 1 and Nf = 2 spectra 
at TTi = 0.15, as shown in flgure ^. And indeed, when we calculate the chiral condensate 



from the Nf = 2 spectral densities, we obtain a condensate, as shown in flgure ^, that is 
quite clearly vanishing as m — » 0. 

We also calculate (Q^), as listed in table |I^. If we plot the values versus m, as in 
figure ^ we see that the dependence is not linear, as one would expect in the phase with 
symmetry breaking, but rather is roughly oc m'^f = m? as one would expect in the chirally 
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symmetric phase. 

Our calculations are thus consistent with the theory being in a phase with explicit chiral 
symmetry. Now, while one does indeed expect to find chiral symmetry restoration as we 



increase Nf, it is not expected to occur before Nj = 6 or so ||20|. So if it already happens 
in our calculations for Nf = 2 this would suggest that our model provides a rather poor 
representation of the fermionic physics in QCD. 

There is another possibility: that for m < 0.3 we have been driven by finite-volume 
effects into a symmetry restored phase. A mass oi m = 0.3 corresponds to m ^ 25MeV 
(using our usual criterion for converting to physical units) and this is surely small enough, 
in a volume V ~ 2.5fm'^, to lead to strong finite size effects. Moreover, such a change of 



phase has been observed in lattice QCD simulations |2T|. To explore this possibility we need 
to perform calculations for larger m. Since we are reluctant to use m ~ Xnz, "we now turn 
to the fixed- A^ calculations. 

4. Nf = 2; fixed N. 

Our fixed N calculations are performed over the range m G [0.15,3.0]. In the mass 
range m < 0.50 we find the same qualitative behaviour as we obtained in the variable-iV 



calculations. However, when we compare the m = 0.15 and m = 3.0 spectra, as in figure |3^, 
we immediately see a great difference between the two: the spectra corresponding to the 
smaller mass shows the depletion of eigenvalues at small A that we expect to see if chiral 



symmetry is to be restored (compare with figure ^9]); the spectrum corresponding to m = 3.0 
shows no such depletion. 

We integrate the spectra for the various masses, and plot the resulting chiral condensate 
in figure ^. This figure appears to confirm our speculation that chiral symmetry restoration 
is a small-m finite volume effect. For rri > 1.0 ~ SOMeV we have a chiral condensate that 
appears to be heading towards a non-zero 171 = intercept. For m < 1.0 there is a sudden and 
drastic suppression. Such a behaviour is what one would expect if the restoration of chiral 
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symmetry were due to finite volume effects. The corresponding plot of {Q'^{m)) is given 



in figure B3. We see an approximately linear behaviour for the larger quark masses, as we 



would expect for chiral symmetry breakdown, and, approximately quadratic behaviour for 



the smaller quark masses as before (see figure pl|). All this strongly suggests that instantons 



do indeed break chiral symmetry for two flavours, but that finite volume effects set in earlier 
than with one flavour. Of course, one really needs to perform such calculations on larger 
volumes, in order to render this plausible scenario completely convincing. 

We note that all the Nj = 2 spectral densities have a power-like divergence at A = 0. 
The parameters of the fits to b/X'^ are listed in table 0. For those values of m where we have 
chiral symmetry breaking the power of the divergence appears to be (i ^ 0.3 as for Nf = 1. 
As we enter the range of m where chiral symmetry is restored, the value of d increases, and 
the coefficient b decreases very rapidly. (Note that this supports our earlier suggestion that 
our Nf = 1 calculation with m = 0.15 is afflicted by similar finite- V^ effects.) There is some 
evidence from the m = 1.0 and 2.0 fits that 6(m)/A'^^™^ = bQ{m/XY° as for Nf = 1. 

Our estimates for the t] mass are listed in table |V|. We see that for values of m where we 
observe the spontaneous breaking of chiral symmetry, the mass is very similar for Nf = 1 
and Nf = 2. The restoration of chiral symmetry is associated with m > becoming much 
heavier (and, indeed, much harder to estimate reliably). 

VI. CONCLUSIONS 

We have used a simple model for the way instanton zero-modes mix with each other, 
in order to learn how topology contributes to the Dirac spectral density, particularly at 
the small eigenvalues that are relevant to chiral symmetry breaking. We have done so for 
quenched QCD and for full QCD with both one and two flavours. 

We have found that instantons do indeed appear to break chiral symmetry spontaneously, 
both in quenched QCD and in full QCD. We also found that the topological susceptibility 
appears to vanish linearly with the quark mass m, as required by the QCD anomalous Ward 
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identities when in a phase with chiral symmetry breaking. It is remarkable that our simple 
model manages to capture the subtle trade-off between exact zero-modes and the chiral 
symmetry breaking modes arbitrarily close to zero, which underlies this result. 

That instantons should contribute to chiral symmetry breaking is not unexpected flj-Q. 
Our second result is less conventional. We find that the Dirac spectral density one obtains 
from instantons has a powerlike divergence as A — » 0. In quenched QCD this ~ l/X^ 
behaviour translates into a pathological l/mf" divergence in the chiral condensate. In the 
case of full QCD the divergence is pushed to smaller eigenvalues as A ^ 0, perhaps having 
the form {m/XY, and the chiral condensate is well-behaved. Nonetheless such a divergence 
means that the A -^ and m ^ limits need special care, thus undermining, for example, 
the usual Banks-Casher relation. 

These conclusions broadly confirm the results of earlier work [0,0] which used (essen- 



tially) the same model. The most important way in which we have extended that work, at 
least for quenched QCD, is to show, by varying the various components of the model, that 
the divergence appears to be a generic feature of instanton mixing. In addition, we have 
also shown that the strength of the divergence depends on the instanton density, becom- 
ing negligible for very high densities. In fact it seems that 'realistic' instanton ensembles, 
as produced via lattice simulations [§], fall precisely in the intermediate range of densities 
where the strength of the divergence is varying most rapidly ||10[|. Thus we are not able 



to predict the precise phenomenological impact of this divergence upon quenched QCD. 
However it is interesting that when we calculate the instanton mixing within these lattice 
ensembles, we find that the smaller instantons contribute disproportionately to the smallest 
eigenvalues; as though they were part of a more dilute gas. This may provide the crucial 
link between the rather dense gases of large instantons that many lattice studies obtain, and 
the much more dilute gases of significantly smaller instantons that appear to be required for 
the phenomenological success of liquid instanton models 0. 

For full QCD the range of our calculations has been necessarily limited and further work 
is needed. We have performed exploratory calculations of some flavour- singlet meson masses 
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and appear to obtain an r( mass that remains massive in the chiral hmit, and is ~ IGeV if 
we introduce 'physical units' in the simplest possible fashion. Mass calculations using the 
eigenvectors of the Dirac operator would be of particular interest, not least because they 
would allow us to test for the presence of a massless Goldstone pion. Moreover, the fact 
that the eigenvectors become much more extended in space-time as A — i> (as we saw in 
quenched QCD) should have phenomenological consequences. 

While we now feel that we have convincing evidence for the claim that the mixing of 
would-be instanton zero-modes will generically produce chiral symmetry breaking and a 
A — > divergence in the Dirac spectral density, we cannot of course be sure that this 
survives in the full field theory. For example, if there is some other dynamical mechanism 
that produces a non-zero density of modes at A = 0, then the latter will mix with the 
(mixed) would-be zero-modes near A = 0, and we cannot be certain that the divergence in 
the latter will survive. Equally, if it should be that the modes in the divergent peak have a 
space-time extent that — > oo as A ^ 0, then they might well be suppressed by the confining 
fluctuations of the theory. There are obvious ways that we can, and should, investigate these 
possibihties within our model. 

An intriguing question is whether there is any relation between the divergence that we 
have found and the logarthmic divergence that appears in quenched chiral perturbation 
theory [|l^. A link is plausible because the strength of the chiral logarithm is determined 
by the strength, 6, of the pseudoscalar flavour singlet annihilation diagram which itself is 
determined by the topological fluctuations of the quenched vacuum (and whose iteration, in 
full QCD, produces the t]' mass). If 6 is indeed small at the large instanton densities where 
our divergence becomes weak enough to be fit with a logarithm [|10|, then the consistency of 
the two approaches becomes possible. Indeed one can argue |T^ that the chiral logarithm 
will become a power if one sums leading logs to all orders in chiral perturbation theory. Of 
course once 6 becomes large these perturbative calculations might become unreliable. In 
that case our model may provide a non-perturbative method for calculating the fate of these 
quenched chiral logarithms. Of course, one might ask how the power divergence that we see 
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in full QCD fits into all this. We conjecture that it might be related, in a fashion analogous 
to the above, to the remnant chiral logarithms in partially quenched QCD. This scenario 
can and should be tested within our model. 
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TABLES 



Type 


/ 


V 


Nc 


d 


xI/Ndf 


XV^DF 


Hard Sphere 


0.2 


1.0 


650000 


0.656±0.003 


2.40 


600 




0.5 


1.0 


128000 


0.695±0.002 


1.92 


508 




1.0 


0.41 


130000 


0.540±0.003 


1.71 


123 




1.0 


1.0 


126000 


0.595±0.002 


1.38 


258 




1.0 


1.52 


97000 


0.617±0.001 


1.23 


400 




1.0 


2.44 


93000 


0.640±0.001 


1.47 


702 




1.0 


7.72 


9780 


0.668±0.003 


1.73 


456 




1.75 


1.0 


111000 


0.309±0.002 


1.94 


48 




2.5 


1.0 


63200 


0.075±0.014 


1.72 


2.59 




5.3 


1.0 


6720 


0.004±0.010 


2.20 


2.19 




10.0 


1.0 


1266 


0.008±0.016 


1.91 


1.90 


Gaussian 


1.0 


1.0 


12600 


0.588±0.005 


1.36 


45 




1.0 


2.4 


15500 


0.634±0.002 


1.31 


201 




2.5 


1.0 


6320 


0.290±0.007 


1.46 


4.87 


Classical 


1.0 


1.0 


2100 


0.477±0.010 


1.32 


5.33 




1.0 


2.4 


450 


0.538±0.022 


1.38 


3.50 


Class. (II) 


1.0 


1.0 


2520 


0.555±0.006 


1.56 


34.6 



TABLE I. Parameters and results for (5 = 0, iVj = ensembles. / is the packing fraction, V 
is the space-time volume, N^ is the number of configurations in the ensemble, d is the degree of 
divergence, Xp/^DF is the chi-square for the best power fit to the divergence and xf/^DF is the 
chi-square for the best log fit. 
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Type 


/ 


V 


Nc 


d 


xI/Ndf 


XV^DF 


Hard Sphere 


1.0 
1.0 


1.0 

2.4 


126000 
62000 


0.363±0.006 
0.522±0.003 


1.90 

2.07 


65 
253 



TABLE II. Parameters and results for Nf = ensembles with variable Q and N. f is now the 
mean packing fraction. 



Type 


/ 


V 


(5/9(0,2) 


M2,4) 


M4,6) 


Hard Sphere 


0.77 
2.96 


0.41 
0.11 


0.88(3) 
0.85(3) 


0.75(4) 
0.75(3) 


0.68(3) 
0.67(3) 



TABLE III. Values of 6p{\Q\,\Q\ + 2), as defined in the text, for two different Nf = hard 
sphere ensembles. 
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Set 


^NZ 


Nf 


m 


{Q') 


{Ni) 


A 


2.0 


1 


0.5 


25.65±0.45 


35.53 


B 


2.0 


1 


0.3 


14.28±0.16 


34.91 


C 


2.0 


1 


0.15 


8.13±0.04 


35.75 


D 


2.0 


2 


0.5 


5.04±0.04 


21.35 


E 


2.0 


2 


0.3 


2.43±0.03 


23.17 


F 


2.0 


2 


0.15 


0.72±0.01 


30.85 


G 


- 




3.0 


98.615±0.78 


62.80 


H 


- 




2.0 


80.097±0.92 


63.03 


I 


- 




1.0 


60.383±0.49 


62.95 


J 


- 




0.5 


34.22±0.26 


62.98 


K 


- 




0.3 


22.00±0.19 


62.99 


L 


- 




0.15 


11.66±0.07 


63.01 


M 


- 


2 


3.0 


83.20±0.64 


63.07 


N 


- 


2 


2.0 


62.49±0.40 


62.90 





- 


2 


1.0 


31.37±0.13 


62.93 


P 


- 


2 


0.5 


ll.61ib0.10 


62.98 


Q 


- 


2 


0.3 


4.71±0.05 


62.96 


R 


- 


2 


0.15 


0.89±0.02 


62.98 



TABLE IV. Some information about the iVj ^ full QCD ensembles. Where N has been kept 
fixed, a value for X^z is not needed, and is not given. 
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Set 


b 


d 


xVNdf 


m^, 


m^ 


A 


2.904±0.158 


0.668±0.008 


1.09 


11.2±2.0 


11.3±2.5 


B 


1.956±0.181 


0.699±0.017 


0.94 


12.9±2.5 


10.9±1.8 


C 


0.491±0.058 


0.861±0.022 


2.41 


11.5±1.5 


10.2±2.4 


D 


0.228±0.031 


0.882±0.021 


1.15 


17.3±1.9 


15.2±2.0 


E 


0.043±0.009 


0.979±0.032 


1.44 


19.3±2.0 


17.4±2.7 


F 


- 


- 


- 


20.0±2.0 


16.7±2.2 


G 


30.917±3.324 


0.300±0.018 


0.99 


8.9±3.1 


- 


H 


28.766±2.058 


0.307±0.013 


1.35 


8.5±2.5 


- 


I 


36.517±4.636 


0.252±0.019 


1.41 


11.5±2.2 


- 


J 


24.048±2.936 


0.288±0.019 


2.05 


11.7±2.1 


- 


K 


17.599±2.632 


0.309±0.023 


3.27 


12.4±2.9 


- 


L 


4.395±0.586 


0.478±0.024 


4.20 


12.1±2.3 


- 


M 


30.410±3.919 


0.296±0.020 


1.17 


9.5±3.7 


- 


N 


30.990±2.520 


0.285±0.014 


1.57 


13.6±2.7 


- 





17.665±1.714 


0.343±0.017 


1.41 


15.0±3.6 


- 


P 


4.404±0.267 


0.519±0.013 


1.99 


20.6±2.5 


- 


Q 


0.580±0.051 


0.788±0.019 


1.55 


15.8±2.5 


- 


R 


0.006±0.002 


1.307±0.043 


2.30 


18.2±2.1 


- 



TABLE V. Parameters of power-law fits to the spectral densities, for the ensembles given in 



table IV, as well as mass estimates for the r] and a mesons. 
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FIG. 1. Exponent of power divergence, d, versus space-time volume V. Hard sphere (o) and 
Gaussian (*) wavef unctions, with density / = 1. 



14 
12 
10 

b 8 
6 
4 
2 



2 3 4 

V 



FIG. 2. Volume dependence of the coefficient b of the 1/A piece in the spectral density. For 
hard spheres as in figure |l| 
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FIG. 3. Exponent of power divergence, d, versus packing fraction /. Hard sphere (o) and 
Gaussian (*) wavef unctions, with volume V = 1. 
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FIG. 4. Spectral density using classical zero mode wavefunctions. Our usual approximation to 
iP[A] (solid) compared to that using the linear addition ansatz (dotted). 
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FIG. 5. The spectral densities obtained at / = 1, V = 1 from hard sphere (sohd), Gaussian 
(o) and classical zero mode (*) wavefunctions. 
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FIG. 6. Comparing spectral densities on R^ (o) and on T^ (solid), with hard sphere wavefunc- 
tions for / = 1, V = 1. 
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FIG. 7. Spectral densities with variable Q on two volumes: V = 1 (solid) and V ~ 2.44 
(dashed). With hard spheres and / = 1. 




FIG. 8. The difference of the two spectral densities in figure 0. 
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FIG. 9. Actual spectral density (solid) compared to background curve (*) from pairwise eigen- 
value splitting as defined in the text. Hard spheres with V = 1, / = 1. 
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FIG. 10. Spectral densities: instanton gas (solid), dipole gas with maximum separation 
Ix"*" — x~| < 2/9 (dotted); background curve for dipole gas (long dashed). Using hard spheres 
and y = 1, / = 1. 
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FIG. 11. Actual spectral density (dotted) compared to background curve (•) for a dipole gas 
with niaxiniuni separation \x~^ — x~\ < 2p. With V = 1, / = 1. 
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FIG. 12. Size of eigenfunctions of iJZ) versus A. For volumes V = 0.4096 {-k), V = 1 (o) and 
V f» 1.5 (solid) using / = 1 hard spheres. 



61 



250 



200 



150 







FIG. 13. The spectral density for / = 1, V = 1, Nf = plotted for the full range of A. 
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FIG. 14. Nf = 1, X^z = 2.0, V = 1, m = 0.15. The total number of objects as a functfon of 
the configuration number in the Monte Carlo sequence. 
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FIG. 15. Nf = 1, Xi\iz = 2.0, V = 1, m = 0.15. The net winding number as a function of the 
configuration number. We show only a section of the run for clarity. 
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FIG. 16. Nf = 1, Xnz = 2.0, V = \, m = 0.5. The total number of objects as a function of 
the configuration number. 
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FIG. 17. Nf = 1, Xpfz = 2.0, V = 1, m = 0.5. The net winding number as a function of the 
configuration number. We show only a section of the run for clarity. 
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FIG. 18. Nf = 1, Xnz = 2.0, V = 1. The spectral density for quark mass m = 0.5 (solid), 0.3 
(long dashed) and 0.15 (dotted). 
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FIG. 19. The quark condensate ('0V') ^^ ^ function of the quark mass m from the spectra 
plotted in figure O. 
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FIG. 20. Nf = 1, \]\fz = 2.0, V = 1. The second moment of the winding number distribution 
(Q^) as a function of the quark mass m. 
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FIG. 21. Nf = 1, Xnz = 2.0, V = 1, m = 0.15. Correlation functions for rj {*), and for a (o). 
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FIG. 22. Nf = 1, Xnz = 2.0, V = 1, m = 0.15. The ry effective mass from the correlation 



functions shown in figure |21 
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FIG. 23. Nf = 1, Fixed N, V = 1. The spectral density for quark mass m = 0.5 (solid), 0.3 
(long dashed) and 0.15 (dotted). 
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FIG. 24. Nf = 1, Fixed N, V = 1. The quark condensate (V'V') &s a function of the quark 



mass m. 
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FIG. 25. Nf = 1, V = 1. The mass of the 77 as a function of the quark mass. 
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FIG. 26. Nf = 2, Aat^ = 2.0, V = 1, m = 0.15. The total number of objects as a function of 
the configuration number. 
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FIG. 27. Nf = 2, Xi\iz = 2.0, F = 1, m = 0.15. The net winding number as a function of the 
configuration number. We show only a section of the run for clarity. 
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FIG. 28. Nf = 2, Xnz = 2.0, V = 1. The spectral density for quark mass m = 0.5 (solid), 0.3 
(long dashed) and 0.15 (dotted). 
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FIG. 29. Aat^ = 2.0, V = 1, m = 0.15. A comparison of the spectral densities obtained for 
Nf = 1 (dark-solid) and Nf = 2 (light-dotted). 
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FIG. 30. The quark condensate (V'V') as a function of the quark mass m from the A^j = 2 
spectra plotted in figure pi. 
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FIG. 31. Nf = 2, X]\fz = 2.0, V = 1. The second moment of the winding number distribution 
(Q^) as a function of the quark mass m. 
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FIG. 32. Nf = 2, Fixed A^, V = 1. The spectral density for quark mass m = 0.15 (sohd) and 
m = 3.0 (long dashed). 
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FIG. 33. Nf = 2, Fixed N, V = 1. The quark condensate (V'V') &s a function of the quark 



mass m. 




FIG. 34. Nf = 2, Fixed N, V = 1. The second moment of the winding number distribution 
(Q^) as a function of the quark mass m. 
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